Подгружаем всё необходимое¶
In [95]:
import json
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import matplotlib.colors as mcolors
import seaborn as sns
from statsmodels.miscmodels.ordinal_model import OrderedModel
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
import re
In [96]:
df = pd.read_csv('DTP_DATA_2025_PROCESSED.csv', index_col=0)
df
Out[96]:
| REGION | DATE | COORD_L | COORD_W | road_name | road_category | n_VEHICLES | n_PARTICIPANTS | ID | n_DEATHS | ... | site_objects_cat | severity | YEAR | MONTH | WEEKDAY | SEASON | is_WEEKEND | HOUR | is_NIGHT | is_PEAK_HOUR | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 6 | 1 | 31.01.2015 | 81.151944 | 53.740000 | Романово - Завьялово - Баево - Камень-на-Оби | 5.0 | 1 | 3 | 161242174 | 0 | ... | 12 | 2 | 2015 | 1 | 5 | 1 | 1 | 9 | 0 | 1 |
| 8 | 1 | 30.01.2015 | 85.018056 | 51.684444 | Куяган - Куяча - Тоурак | 6.0 | 2 | 3 | 161105683 | 0 | ... | 12 | 2 | 2015 | 1 | 4 | 1 | 0 | 14 | 0 | 0 |
| 12 | 1 | 30.01.2015 | 81.250000 | 53.818056 | Барнаул - Камень-на-Оби - граница Новосибирско... | 5.0 | 2 | 3 | 161763431 | 0 | ... | 12 | 1 | 2015 | 1 | 4 | 1 | 0 | 17 | 0 | 1 |
| 39 | 1 | 24.01.2015 | 51.000000 | 84.000000 | Быканов Мост - Солоновка - Солонешное - границ... | 7.0 | 1 | 2 | 160331994 | 0 | ... | 12 | 1 | 2015 | 1 | 5 | 1 | 1 | 19 | 0 | 0 |
| 42 | 1 | 23.01.2015 | 84.000000 | 53.000000 | Быканов Мост - Солоновка - Солонешное - границ... | 7.0 | 1 | 2 | 160213415 | 1 | ... | 12 | 3 | 2015 | 1 | 4 | 1 | 0 | 21 | 0 | 0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 1475642 | 10011 | 18.12.2022 | 55.731164 | 67.417464 | г. Нарьян-Мар - г. Усинск, участок г. Нарьян-М... | 5.0 | 1 | 1 | 222679998 | 1 | ... | 3 | 3 | 2022 | 12 | 6 | 1 | 1 | 15 | 0 | 0 |
| 1475665 | 10011 | 16.06.2024 | 54.095821 | 67.629532 | г. Нарьян-Мар - г. Усинск, участок г. Нарьян-М... | 6.0 | 1 | 2 | 223892751 | 0 | ... | 12 | 1 | 2024 | 6 | 6 | 3 | 1 | 9 | 0 | 1 |
| 1475669 | 10011 | 04.07.2024 | 53.607541 | 67.836463 | Подъездная автомобильная дорога к п. Красное | 1.0 | 1 | 1 | 223948754 | 0 | ... | 6 | 2 | 2024 | 7 | 3 | 3 | 0 | 9 | 0 | 1 |
| 1475671 | 10011 | 25.08.2024 | 53.612724 | 67.837285 | Подъездная автомобильная дорога к п. Красное | 4.0 | 1 | 1 | 224161731 | 0 | ... | 12 | 1 | 2024 | 8 | 6 | 3 | 1 | 14 | 0 | 0 |
| 1475672 | 10011 | 24.09.2024 | 55.956877 | 67.303057 | г. Нарьян-Мар - г. Усинск, участок п. Харьягин... | 8.0 | 2 | 2 | 224282217 | 0 | ... | 12 | 2 | 2024 | 9 | 1 | 4 | 0 | 20 | 0 | 0 |
473188 rows × 79 columns
Предобработка¶
In [97]:
df = df.drop(columns=['n_pedestrians', 'russian_vehicle'])
In [98]:
df = df.rename(columns={'vehicle_age_max': 'vehicle_year_max', 'vehicle_age_min': 'vehicle_year_min', 'vehicle_age_avg': 'vehicle_year_avg'})
Удобный формат регионов¶
In [99]:
regions = pd.read_json('regions.json')
regions
Out[99]:
| id | name | districts | |
|---|---|---|---|
| 0 | 71140 | Ямало-Ненецкий АО | [{"id": "71166", "name": "Шурышкарский район"}... |
| 1 | 71100 | Ханты-Мансийский АО | [{"id": "71126", "name": "Сургутский район"}, ... |
| 2 | 71 | Тюменская область | [{"id": "712581", "name": "Ярковский район"}, ... |
| 3 | 10011 | Ненецкий АО | [{"id": "11100", "name": "Ненецкий АО"}] |
| 4 | 99 | Еврейская автономная область | [{"id": "99205", "name": "Биробиджанский район... |
| ... | ... | ... | ... |
| 80 | 4 | Красноярский край | [{"id": "4233", "name": "Минусинский район"}, ... |
| 81 | 3 | Краснодарский край | [{"id": "3405", "name": "г. Армавир"}, {"id": ... |
| 82 | 1 | Алтайский край | [{"id": "14011", "name": "г.Барнаул"}, {"id": ... |
| 83 | 35 | Республика Крым | [{"id": "8125", "name": "Ялта"}, {"id": "8124"... |
| 84 | 67 | г. Севастополь | [{"id": "8126", "name": "Севастополь"}, {"id":... |
85 rows × 3 columns
In [100]:
regions['id'] = regions['id'].astype(str)
region_map = dict(zip(regions['id'], regions['name']))
df['region_name'] = df['REGION'].astype(str).map(region_map)
In [101]:
df['region_name']
Out[101]:
6 Алтайский край
8 Алтайский край
12 Алтайский край
39 Алтайский край
42 Алтайский край
...
1475642 Ненецкий АО
1475665 Ненецкий АО
1475669 Ненецкий АО
1475671 Ненецкий АО
1475672 Ненецкий АО
Name: region_name, Length: 473188, dtype: object
In [102]:
df.columns
Out[102]:
Index(['REGION', 'DATE', 'COORD_L', 'COORD_W', 'road_name', 'road_category',
'n_VEHICLES', 'n_PARTICIPANTS', 'ID', 'n_DEATHS', 'n_INJURED',
'vehicle_failure', 'non_private_vehicle', 'white_vehicle',
'black_vehicle', 'colored_vehicle', 'drunk_driver', 'female_driver',
'escaped', 'no_seatbelt_injury', 'n_drunk', 'n_children', 'n_cyclists',
'vehicle_year_min', 'vehicle_year_max', 'vehicle_year_avg', 'n_class_a',
'n_class_b', 'n_class_c', 'n_class_d', 'n_class_e', 'n_class_s',
'n_front_drive', 'n_rear_drive', 'n_4wd', 'n_guilty', 'guilty_share',
'n_fatal_violations', 'guilty_exp_avg', 'exp_avg', 'road_rank_cat',
'road_defects_cat', 'traffic_changes_bin', 'traffic_changes_cat',
'road_surface_cat', 'TYPE_cat', 'out_of_town', 'street_rank_cat',
'weather_interpretable', 'weather_cat', 'adj_objects_interpretable',
'adj_objects_cat', 'cause_factors_cat', 'crossing_violation',
'impaired_driving', 'interference_violation', 'license_violation',
'maneuver_violation', 'other_violation', 'pedestrian_violation',
'sudden_appearance_violation', 'traffic_control_violation',
'transport_violation', 'vehicle_tech_violation', 'wrong_way',
'no_lighting', 'lighting_cat', 'site_objects_cat', 'severity', 'YEAR',
'MONTH', 'WEEKDAY', 'SEASON', 'is_WEEKEND', 'HOUR', 'is_NIGHT',
'is_PEAK_HOUR', 'region_name'],
dtype='object')
Освещение¶
In [103]:
light_map = {0: 'День', 1: 'Не указано', 2: 'Ночь без освещения', 3: 'Ночь с освещением', 4: 'Сумерки'}
df['lighting'] = df['lighting_cat'].astype(int).map(light_map)
df['lighting']
Out[103]:
6 Сумерки
8 День
12 День
39 Ночь без освещения
42 Ночь без освещения
...
1475642 Ночь без освещения
1475665 День
1475669 День
1475671 День
1475672 Ночь без освещения
Name: lighting, Length: 473188, dtype: object
Группируем регионы в округа¶
In [104]:
with open('regions_federal_districts_89.json', 'r', encoding='utf-8') as f:
region_map = json.load(f)
df['district'] = df['region_name'].map(region_map)
Графики до очистки¶
In [105]:
try:
plt.rcParams['font.family'] = 'Montserrat'
except:
plt.rcParams['font.family'] = 'sans-serif'
sns.set_style("white")
# 2. Данные
# ВАЖНО: Считаем квартили по ВСЕМ данным (до фильтрации),
# чтобы выбросы не исказили статистику, или по уже очищенным - зависит от вашей логики.
# Обычно IQR считается на исходном массиве.
raw_data = df['guilty_exp_avg'].dropna()
# --- РАСЧЕТ ГРАНИЦЫ (IQR METHOD) ---
Q1 = raw_data.quantile(0.25)
Q3 = raw_data.quantile(0.75)
IQR = Q3 - Q1
cutoff_value = Q3 + 1.5 * IQR
# -----------------------------------
# 3. Настройка фигуры
plt.figure(figsize=(14, 8), dpi=150)
# 4. Строим график
ax = sns.histplot(
data=raw_data,
bins=60, # Чуть больше бинов, чтобы видеть детали
kde=True,
color="#0261C7", # Спокойный сине-голубой (как на скриншоте)
alpha=0.6,
edgecolor="white",
line_kws={'linewidth': 2.5, 'color': '#153AE0'} # Синяя линия тренда
)
# 5. Линии
# Медиана (Оранжевая)
median_val = raw_data.median()
plt.axvline(median_val, color='#B50E3B', linestyle='--', linewidth=2.5,
label=f'Медиана: {median_val:.1f} лет')
# --- ЛИНИЯ ОТСЕЧЕНИЯ (IQR) ---
plt.axvline(cutoff_value, color='#2d3436', linestyle=':', linewidth=3,
label=f'Граница выбросов: {cutoff_value:.1f} лет')
# (Опционально) Заштриховать зону выбросов
plt.axvspan(cutoff_value, raw_data.max(), color='#B50E3B', alpha=0.05, label='Зона отсечения')
# -----------------------------
# 6. Форматтер оси Y (тысячи)
def y_fmt(x, pos):
return f'{x/1000:.0f}'
ax.yaxis.set_major_formatter(ticker.FuncFormatter(y_fmt))
# 7. Декор
sns.despine(left=False, bottom=False) # Или как тебе удобнее
# Сетка пунктирная
ax.yaxis.grid(True, linestyle='--', which='major', color='grey', alpha=0.3)
ax.set_axisbelow(True)
# Подписи
plt.title('Распределение водительского стажа виновников', fontsize=24, fontweight='bold', color='#2d3436', pad=25)
plt.xlabel('Стаж вождения (лет)', fontsize=18, color='#636e72', labelpad=14)
plt.ylabel('Количество (тыс.)', fontsize=18, color='#636e72', labelpad=14)
# Легенда
plt.legend(frameon=False, fontsize=16, loc='upper right')
# Ограничим X, чтобы график не улетал в бесконечность, если есть ошибки типа "1900 лет стажа"
# Но покажем хвост с выбросами, чтобы было видно, ЧТО мы отрезаем
plt.xlim(0, max(cutoff_value * 1.5, 110))
plt.tight_layout()
plt.show()
Очистка¶
In [107]:
exp = df['exp_avg']
exp_IQR = exp.quantile(0.75) - exp.quantile(0.25)
exp_edge = exp.quantile(0.75) + 1.5 * exp_IQR
exp_edge
Out[107]:
48.0
In [108]:
guilty_exp = df['guilty_exp_avg']
guilty_exp_IQR = guilty_exp.quantile(0.75) - guilty_exp.quantile(0.25)
guilty_exp_edge = guilty_exp.quantile(0.75) + 1.5 * guilty_exp_IQR
guilty_exp_edge
Out[108]:
56.0
In [109]:
df = df[(df['guilty_exp_avg'] <= guilty_exp_edge) & (df['exp_avg'] < guilty_exp_edge)]
df
Out[109]:
| REGION | DATE | COORD_L | COORD_W | road_name | road_category | n_VEHICLES | n_PARTICIPANTS | ID | n_DEATHS | ... | MONTH | WEEKDAY | SEASON | is_WEEKEND | HOUR | is_NIGHT | is_PEAK_HOUR | region_name | lighting | district | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 6 | 1 | 31.01.2015 | 81.151944 | 53.740000 | Романово - Завьялово - Баево - Камень-на-Оби | 5.0 | 1 | 3 | 161242174 | 0 | ... | 1 | 5 | 1 | 1 | 9 | 0 | 1 | Алтайский край | Сумерки | Сибирский |
| 8 | 1 | 30.01.2015 | 85.018056 | 51.684444 | Куяган - Куяча - Тоурак | 6.0 | 2 | 3 | 161105683 | 0 | ... | 1 | 4 | 1 | 0 | 14 | 0 | 0 | Алтайский край | День | Сибирский |
| 12 | 1 | 30.01.2015 | 81.250000 | 53.818056 | Барнаул - Камень-на-Оби - граница Новосибирско... | 5.0 | 2 | 3 | 161763431 | 0 | ... | 1 | 4 | 1 | 0 | 17 | 0 | 1 | Алтайский край | День | Сибирский |
| 39 | 1 | 24.01.2015 | 51.000000 | 84.000000 | Быканов Мост - Солоновка - Солонешное - границ... | 7.0 | 1 | 2 | 160331994 | 0 | ... | 1 | 5 | 1 | 1 | 19 | 0 | 0 | Алтайский край | Ночь без освещения | Сибирский |
| 42 | 1 | 23.01.2015 | 84.000000 | 53.000000 | Быканов Мост - Солоновка - Солонешное - границ... | 7.0 | 1 | 2 | 160213415 | 1 | ... | 1 | 4 | 1 | 0 | 21 | 0 | 0 | Алтайский край | Ночь без освещения | Сибирский |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 1475635 | 10011 | 06.08.2022 | 53.104005 | 67.671952 | Автомобильная дорога по улице Угольная | 6.0 | 1 | 2 | 222381293 | 0 | ... | 8 | 5 | 3 | 1 | 5 | 1 | 0 | Ненецкий АО | День | Северо-Западный |
| 1475642 | 10011 | 18.12.2022 | 55.731164 | 67.417464 | г. Нарьян-Мар - г. Усинск, участок г. Нарьян-М... | 5.0 | 1 | 1 | 222679998 | 1 | ... | 12 | 6 | 1 | 1 | 15 | 0 | 0 | Ненецкий АО | Ночь без освещения | Северо-Западный |
| 1475665 | 10011 | 16.06.2024 | 54.095821 | 67.629532 | г. Нарьян-Мар - г. Усинск, участок г. Нарьян-М... | 6.0 | 1 | 2 | 223892751 | 0 | ... | 6 | 6 | 3 | 1 | 9 | 0 | 1 | Ненецкий АО | День | Северо-Западный |
| 1475669 | 10011 | 04.07.2024 | 53.607541 | 67.836463 | Подъездная автомобильная дорога к п. Красное | 1.0 | 1 | 1 | 223948754 | 0 | ... | 7 | 3 | 3 | 0 | 9 | 0 | 1 | Ненецкий АО | День | Северо-Западный |
| 1475672 | 10011 | 24.09.2024 | 55.956877 | 67.303057 | г. Нарьян-Мар - г. Усинск, участок п. Харьягин... | 8.0 | 2 | 2 | 224282217 | 0 | ... | 9 | 1 | 4 | 0 | 20 | 0 | 0 | Ненецкий АО | Ночь без освещения | Северо-Западный |
419595 rows × 80 columns
In [110]:
len(df)
Out[110]:
419595
После очистки осталось 419595 строк¶
Графички после очистки¶
In [111]:
hex_colors = ['#B50E3B', '#850AD6', '#490FD2', '#153AE0', '#0261C7']
# 2. Подготовка данных
# Сортируем от большего к меньшему
region_counts = df['region_name'].value_counts()
n_bars = len(region_counts)
# 3. Создаем градиентную палитру на N регионов
custom_cmap = mcolors.LinearSegmentedColormap.from_list("custom_gradient", hex_colors, N=n_bars)
gradient_palette = custom_cmap(np.linspace(0, 1, n_bars))
# 4. Настройка шрифта
try:
plt.rcParams['font.family'] = 'Montserrat'
except:
plt.rcParams['font.family'] = 'sans-serif'
# 5. Рисуем график
# ВАЖНО: figsize=(15, 30) — делаем график очень длинным по вертикали!
plt.figure(figsize=(15, 30), dpi=150)
sns.set_style("white")
ax = sns.barplot(
x=region_counts.values,
y=region_counts.index,
hue=region_counts.index,
palette=gradient_palette,
edgecolor='none',
legend=False
)
# 6. Оформление
sns.despine(left=True, bottom=False)
plt.title('Количество ДТП по всем регионам РФ', fontsize=24, fontweight='bold', color='#2d3436', pad=30, loc='left')
plt.xlabel('Количество ДТП (тыс.)', fontsize=14, color='gray', labelpad=15)
plt.ylabel('', fontsize=12)
# Увеличим шрифт названий регионов, чтобы читалось
ax.tick_params(axis='y', labelsize=11)
# 7. Форматирование оси X (тысячи)
ax.xaxis.set_major_formatter(ticker.FuncFormatter(lambda x, pos: f'{x/1000:.0f}'))
ax.xaxis.grid(True, linestyle='--', color='grey', alpha=0.2)
ax.xaxis.tick_top() # Перенесем шкалу наверх, чтобы было удобнее смотреть длинный список
ax.xaxis.set_label_position('top')
# 8. Цифры (только если столбец достаточно большой, чтобы не мусорить)
for i, v in enumerate(region_counts.values):
# Добавляем подпись
ax.text(
v + (region_counts.max() * 0.01), # Небольшой отступ
i,
f'{v:,.0f}'.replace(',', ' '),
va='center',
fontsize=10,
color='#636e72'
)
plt.tight_layout()
plt.show()
/var/folders/5y/prlkw2jx437djr2ybr7mqrfm0000gn/T/ipykernel_12109/3431987402.py:23: UserWarning: Numpy array is not a supported type for `palette`. Please convert your palette to a list. This will become an error in v0.14 ax = sns.barplot(
In [112]:
region_counts = df['region_name'].value_counts()
top_5 = region_counts.head(5)
bottom_5 = region_counts.tail(5)
# Единый масштаб
global_max = top_5.max() * 1.15
# 2. Палитра
hex_colors = ['#B50E3B', '#850AD6', '#490FD2', '#153AE0', '#0261C7']
full_cmap = mcolors.LinearSegmentedColormap.from_list("custom", hex_colors, N=100)
colors_top = full_cmap(np.linspace(0, 0.45, 5))
colors_bottom = full_cmap(np.linspace(0.55, 1, 5))
# 3. Настройки шрифтов
plt.rcParams['font.family'] = 'sans-serif'
plt.rcParams['font.sans-serif'] = ['Arial', 'Helvetica', 'Verdana']
# 4. Фигура
fig, axes = plt.subplots(2, 1, figsize=(16, 12), dpi=150)
plt.subplots_adjust(hspace=0.35)
# --- ГРАФИКИ ---
sns.barplot(x=top_5.values, y=top_5.index, hue=top_5.index, palette=colors_top, ax=axes[0], legend=False)
axes[0].set_title('ТОП-5 Регионов по количеству ДТП', fontsize=24, fontweight='heavy', color='#333333', pad=15)
sns.barplot(x=bottom_5.values, y=bottom_5.index, hue=bottom_5.index, palette=colors_bottom, ax=axes[1], legend=False)
axes[1].set_title('5 Регионов с минимальным количеством ДТП', fontsize=24, fontweight='heavy', color='#333333', pad=15)
# --- СТИЛИЗАЦИЯ ---
def style_ax(ax):
ax.set_xlim(0, global_max)
sns.despine(ax=ax, left=True, bottom=False)
ax.set_ylabel('')
# Увеличил шрифт подписи оси
ax.set_xlabel('Количество записей (в тысячах)', fontsize=20, color='gray', labelpad=15)
# Форматтер (тысячи 'k')
ax.xaxis.set_major_formatter(ticker.FuncFormatter(lambda x, pos: f'{x/1000:.0f}'))
# === ИЗМЕНЕНИЕ ЗДЕСЬ: Размер шрифта шкалы (labelsize) ===
ax.tick_params(axis='x', labelsize=18, labelcolor='gray')
ax.xaxis.grid(True, linestyle='--', color='grey', alpha=0.3)
# Названия регионов
ax.tick_params(axis='y', length=0)
plt.setp(ax.get_yticklabels(), fontsize=18, fontweight='bold', color='black')
style_ax(axes[0])
style_ax(axes[1])
# --- ЦИФРЫ У СТОЛБЦОВ ---
def add_labels(ax, data):
for i, v in enumerate(data.values):
ax.text(
v + (global_max * 0.01),
i,
f'{v:,.0f}'.replace(',', ' '),
va='center',
fontsize=18,
fontweight='bold',
color='black'
)
add_labels(axes[0], top_5)
add_labels(axes[1], bottom_5)
plt.show()
/var/folders/5y/prlkw2jx437djr2ybr7mqrfm0000gn/T/ipykernel_12109/2552072259.py:23: UserWarning: Numpy array is not a supported type for `palette`. Please convert your palette to a list. This will become an error in v0.14 sns.barplot(x=top_5.values, y=top_5.index, hue=top_5.index, palette=colors_top, ax=axes[0], legend=False) /var/folders/5y/prlkw2jx437djr2ybr7mqrfm0000gn/T/ipykernel_12109/2552072259.py:26: UserWarning: Numpy array is not a supported type for `palette`. Please convert your palette to a list. This will become an error in v0.14 sns.barplot(x=bottom_5.values, y=bottom_5.index, hue=bottom_5.index, palette=colors_bottom, ax=axes[1], legend=False)
In [113]:
pop_df = pd.read_csv('regions_people.csv', index_col=0)
accidents = df['region_name'].value_counts().reset_index()
accidents.columns = ['region_name', 'accident_count']
# Объединяем таблицы (LEFT JOIN, чтобы не потерять регионы ДТП, но лучше INNER если уверены в названиях)
# В CSV колонка 'name', в df 'region_name'
merged = pd.merge(accidents, pop_df, left_on='region_name', right_on='name', how='inner')
# Считаем относительную величину: (ДТП / Население) * 1000
merged['rate_per_1000'] = (merged['accident_count'] / merged['man']) * 1000
# Сортируем
merged_sorted = merged.sort_values(by='rate_per_1000', ascending=False)
top_5 = merged_sorted.head(5)
bottom_5 = merged_sorted.tail(5)
# --- 3. ВИЗУАЛИЗАЦИЯ (Финальный стиль) ---
# Настройки
plt.rcParams['font.family'] = 'sans-serif'
plt.rcParams['font.sans-serif'] = ['Arial', 'Helvetica', 'Verdana']
hex_colors = ['#B50E3B', '#850AD6', '#490FD2', '#153AE0', '#0261C7']
full_cmap = mcolors.LinearSegmentedColormap.from_list("custom", hex_colors, N=100)
colors_top = full_cmap(np.linspace(0, 0.45, 5))
colors_bottom = full_cmap(np.linspace(0.55, 1, 5))
# Фигура
fig, axes = plt.subplots(2, 1, figsize=(16, 12), dpi=150)
plt.subplots_adjust(hspace=0.4) # Чуть больше места между заголовком и графиком
# Определяем масштаб (здесь лучше не делать единый, если разрыв огромен, но для честности оставим)
# Если в топе цифра 5.0, а внизу 0.2, то единый масштаб "убьет" нижний график.
# Давай сделаем адаптивный масштаб, но визуально схожий стиль.
# ИЛИ, если хочешь показать разрыв, берем max от топа.
global_max = top_5['rate_per_1000'].max() * 1.2
# --- ГРАФИКИ ---
sns.barplot(data=top_5, x='rate_per_1000', y='region_name', palette=colors_top, ax=axes[0], hue='region_name', legend=False)
axes[0].set_title('ТОП-5 Самых опасных регионов (на 1000 жителей)', fontsize=24, fontweight='heavy', color='#333333', pad=15)
sns.barplot(data=bottom_5, x='rate_per_1000', y='region_name', palette=colors_bottom, ax=axes[1], hue='region_name', legend=False)
axes[1].set_title('5 Самых безопасных регионов (на 1000 жителей)', fontsize=24, fontweight='heavy', color='#333333', pad=15)
# --- ОФОРМЛЕНИЕ ---
def style_ax(ax):
# Оставляем единый масштаб для честности сравнения (можно убрать, если низ совсем плоский)
ax.set_xlim(0, global_max)
sns.despine(ax=ax, left=True, bottom=False)
ax.set_ylabel('')
ax.set_xlabel('ДТП на 1 000 человек', fontsize=18, color='gray', labelpad=15)
# Сетка и шрифт шкалы
ax.xaxis.grid(True, linestyle='--', color='grey', alpha=0.3)
ax.tick_params(axis='x', labelsize=16, labelcolor='gray')
# Названия регионов
ax.tick_params(axis='y', length=0)
plt.setp(ax.get_yticklabels(), fontsize=18, fontweight='bold', color='black')
style_ax(axes[0])
style_ax(axes[1])
# --- ЦИФРЫ ---
def add_labels(ax, subset):
for i, row in enumerate(subset.itertuples()):
val = row.rate_per_1000
ax.text(
val + (global_max * 0.01),
i,
f'{val:.2f}', # 2 знака после запятой
va='center',
fontsize=18,
fontweight='bold',
color='black'
)
add_labels(axes[0], top_5)
add_labels(axes[1], bottom_5)
plt.show()
/var/folders/5y/prlkw2jx437djr2ybr7mqrfm0000gn/T/ipykernel_12109/1205570752.py:40: UserWarning: Numpy array is not a supported type for `palette`. Please convert your palette to a list. This will become an error in v0.14 sns.barplot(data=top_5, x='rate_per_1000', y='region_name', palette=colors_top, ax=axes[0], hue='region_name', legend=False) /var/folders/5y/prlkw2jx437djr2ybr7mqrfm0000gn/T/ipykernel_12109/1205570752.py:43: UserWarning: Numpy array is not a supported type for `palette`. Please convert your palette to a list. This will become an error in v0.14 sns.barplot(data=bottom_5, x='rate_per_1000', y='region_name', palette=colors_bottom, ax=axes[1], hue='region_name', legend=False)
In [114]:
plt.rcParams['font.family'] = 'sans-serif'
plt.rcParams['font.sans-serif'] = ['Arial', 'Helvetica', 'Verdana']
sns.set_style("white")
# 2. Подготовка данных
# Сортируем округа по общему количеству ДТП (чтобы слева были самые крупные)
district_order = df['district'].value_counts().index
# 3. Палитра
# Нам нужно 3 цвета (для Severity 1, 2, 3) из твоей градиентной палитры.
# Твоя палитра идет от Малинового (Danger) к Синему (Safe).
# Логично сделать: 3 (Смерть) = Малиновый, 1 (Легкие) = Синий.
hex_colors = ['#B50E3B', '#850AD6', '#490FD2', '#153AE0', '#0261C7']
full_cmap = mcolors.LinearSegmentedColormap.from_list("custom", hex_colors, N=100)
# Берем 3 цвета: Начало (красный), Середина (фиолетовый), Конец (синий)
# Но так как hue обычно сортируется 1->2->3, а мы хотим 3=Красный, нам нужно перевернуть порядок цветов
# или просто взять их вручную.
# Сделаем так: 1 (Легкие) - Синий, 2 (Травмы) - Фиолетовый, 3 (Тяжкие) - Малиновый
severity_palette = [hex_colors[-1], hex_colors[2], hex_colors[0]] # Синий, Фиолетовый, Красный
# 4. Фигура
plt.figure(figsize=(16, 10), dpi=150)
# 5. Строим график
# Используем countplot, он сам посчитает количество
ax = sns.countplot(
data=df,
x='district',
hue='severity',
order=district_order,
palette=severity_palette,
edgecolor='none' # Убираем обводку для чистоты
)
# 6. Оформление
# Заголовок
plt.title('Структура тяжести ДТП по Федеральным Округам', fontsize=24, fontweight='heavy', color='#333333', pad=25)
# Оси
plt.xlabel('Федеральный округ', fontsize=18, fontweight='bold', color='gray', labelpad=15)
plt.ylabel('Количество ДТП (в тыс.)', fontsize=18, fontweight='bold', color='gray', labelpad=15)
# Ось Y (Тысячи)
ax.yaxis.set_major_formatter(ticker.FuncFormatter(lambda x, pos: f'{x/1000:.0f}'))
ax.yaxis.grid(True, linestyle='--', color='grey', alpha=0.3)
ax.tick_params(axis='y', labelsize=14, labelcolor='black')
# Ось X (Округа)
# Поворачиваем текст, так как названия округов длинные
plt.xticks(rotation=25, ha='right', fontsize=14, fontweight='bold', color='black')
# Убираем рамки
sns.despine(left=True, bottom=False)
# 7. Легенда
# Делаем её красивой и понятной
leg = plt.legend(
title='Тяжесть ДТП',
title_fontsize=14,
fontsize=12,
frameon=False,
loc='upper right'
)
# Меняем названия в легенде (если у тебя severity 1, 2, 3)
new_labels = ['1: Легкие', '2: С пострадавшими', '3: Тяжкие/Смерть']
for t, l in zip(leg.texts, new_labels):
t.set_text(l)
plt.tight_layout()
plt.show()
Математическая модель¶
Подготовка датасета¶
In [21]:
df['road_category'] = df['road_category'].astype(int)
In [22]:
# 1. Списки признаков
categorical_cols = [
'district', 'road_defects_cat',
'road_surface_cat', 'lighting_cat', 'SEASON'
]
numerical_cols = [
'n_VEHICLES', 'n_guilty', 'guilty_exp_avg',
'vehicle_failure', 'female_driver', 'no_seatbelt_injury',
'impaired_driving', 'wrong_way'
# Эти столбцы обычно бинарные (0/1), их можно считать численными для модели
]
target = 'severity'
# 2. Собираем X и y
# Берем только нужные колонки и сразу удаляем пустые строки
model_data = df[categorical_cols + numerical_cols + [target]].dropna()
y = model_data[target]
X_cat = model_data[categorical_cols].astype(str) # Категории в строки
X_num = model_data[numerical_cols].astype(float) # Числа во флоат
# 3. One-Hot Encoding для категорий
# drop_first=False, чтобы мы могли сами выбрать базы вручную
X_dummies = pd.get_dummies(X_cat, prefix_sep='_', drop_first=False)
# 4. Объединяем числа и dummy-категории
X = pd.concat([X_num, X_dummies], axis=1)
print(X.columns)
# 5. --- ВЫБОР БАЗОВЫХ КАТЕГОРИЙ ---
# Удаляем столбцы, с которыми будем сравнивать.
# Для новых полей выберем логичные базы (например, самый частый округ или Лето)
cols_to_drop = [
# Старые базы
'lighting_cat_0', # Светлое время
'road_surface_cat_7', # Сухое покрытие
'road_defects_cat_5', # Недостатки содержания (или Нет дефектов)
'SEASON_3', # Сравниваем Зиму/Осень с Летом
'district_Центральный',
]
# Фильтруем, удаляем только те, что реально есть в данных
existing_cols_to_drop = [c for c in cols_to_drop if c in X.columns]
X = X.drop(columns=existing_cols_to_drop)
print(f"Итоговая размерность X: {X.shape}")
print(f"Базовые категории (удалены): {existing_cols_to_drop}")
Index(['n_VEHICLES', 'n_guilty', 'guilty_exp_avg', 'vehicle_failure',
'female_driver', 'no_seatbelt_injury', 'impaired_driving', 'wrong_way',
'district_Дальневосточный', 'district_Приволжский',
'district_Северо-Западный', 'district_Северо-Кавказский',
'district_Сибирский', 'district_Уральский', 'district_Центральный',
'district_Южный', 'road_defects_cat_0', 'road_defects_cat_1',
'road_defects_cat_10', 'road_defects_cat_11', 'road_defects_cat_12',
'road_defects_cat_13', 'road_defects_cat_14', 'road_defects_cat_2',
'road_defects_cat_3', 'road_defects_cat_4', 'road_defects_cat_5',
'road_defects_cat_6', 'road_defects_cat_7', 'road_defects_cat_8',
'road_defects_cat_9', 'road_surface_cat_0', 'road_surface_cat_1',
'road_surface_cat_2', 'road_surface_cat_3', 'road_surface_cat_4',
'road_surface_cat_5', 'road_surface_cat_6', 'road_surface_cat_7',
'lighting_cat_0', 'lighting_cat_1', 'lighting_cat_2', 'lighting_cat_3',
'lighting_cat_4', 'SEASON_1', 'SEASON_2', 'SEASON_3', 'SEASON_4'],
dtype='object')
Итоговая размерность X: (419576, 43)
Базовые категории (удалены): ['lighting_cat_0', 'road_surface_cat_7', 'road_defects_cat_5', 'SEASON_3', 'district_Центральный']
In [23]:
scaler = StandardScaler()
# Масштабируем только численные колонки
num_cols_to_scale = ['n_VEHICLES', 'n_guilty', 'guilty_exp_avg'] # добавь сюда свои численные
# Перезаписываем их в X
X[num_cols_to_scale] = scaler.fit_transform(X[num_cols_to_scale])
Проверка требований¶
Используем порядковую линейную регрессию, target-переменную берём severity, проверяем выполнение ограничения - пропорциональные шансы
In [26]:
try:
plt.rcParams['font.family'] = 'Montserrat'
except:
plt.rcParams['font.family'] = 'sans-serif'
plt.rcParams['font.sans-serif'] = ['Arial', 'Helvetica', 'Verdana']
sns.set_style("white")
y_binary_1 = y.apply(lambda x: 0 if x == 1 else 1)
# Модель Б: Граница между {1,2} и 3
# "Отличает ли фактор Выживание (Легкие + Травмы) от Смерти?"
# 0 = Severity 1 или 2
# 1 = Severity 3
y_binary_2 = y.apply(lambda x: 0 if x <= 2 else 1)
# --- 2. Обучение двух обычных логистических регрессий ---
# Используем тот же X, что и для основной модели!
# class_weight='balanced' обязателен, так как смертельных случаев мало
print("Обучаем бинарную модель 1 (Легкие vs Травмы/Смерть)...")
model_1 = LogisticRegression(max_iter=5000, solver='lbfgs', class_weight='balanced')
model_1.fit(X, y_binary_1)
print("Обучаем бинарную модель 2 (Выживание vs Смерть)...")
model_2 = LogisticRegression(max_iter=5000, solver='lbfgs', class_weight='balanced')
model_2.fit(X, y_binary_2)
# --- 3. Сравнение коэффициентов ---
# Собираем все в один DataFrame
comparison_df = pd.DataFrame({
'Коэфф (1 vs 2+3)': model_1.coef_[0],
'Коэфф (1+2 vs 3)': model_2.coef_[0]
}, index=X.columns)
# Считаем разницу (по модулю)
comparison_df['Разница'] = (comparison_df['Коэфф (1 vs 2+3)'] - comparison_df['Коэфф (1+2 vs 3)']).abs()
# Сортируем: сверху те, кто сильнее всего нарушает правило
comparison_df = comparison_df.sort_values(by='Разница', ascending=False)
# 2. Данные (Те же, что и раньше)
threshold = 0.5
outliers = comparison_df[comparison_df['Разница'] > threshold].copy()
normals = comparison_df[comparison_df['Разница'] <= threshold].copy()
# 3. Фигура (Чуть больше, чтобы влез крупный текст)
plt.figure(figsize=(12, 12), dpi=150)
# 4. Диагональ
lims = [
-0.75,
np.max([comparison_df.max().max(), 1])
]
plt.plot(lims, lims, color='#b2bec3', linestyle='--', linewidth=2.5, label='Идеальная пропорциональность', zorder=1)
# 5. Точки
sns.scatterplot(
x=normals['Коэфф (1 vs 2+3)'],
y=normals['Коэфф (1+2 vs 3)'],
color='#0261C7', alpha=0.5, s=100, edgecolor='white', linewidth=0.5, zorder=2
)
sns.scatterplot(
x=outliers['Коэфф (1 vs 2+3)'],
y=outliers['Коэфф (1+2 vs 3)'],
color='#B50E3B', alpha=1, s=150, edgecolor='white', linewidth=1, zorder=3, label='Специфичные факторы'
)
# 6. Подписи точек (КРУПНЕЕ)
for idx, row in outliers.iterrows():
x_pos = row['Коэфф (1 vs 2+3)']
y_pos = row['Коэфф (1+2 vs 3)']
offset_x = 0.03 if x_pos < y_pos else 0.06
offset_y = 0.03 if x_pos < y_pos else -0.06
plt.text(
x_pos + offset_x,
y_pos + offset_y,
idx,
fontsize=15, # <--- Было 10, стало 13
fontweight='bold',
color='#B50E3B',
zorder=4
)
# 7. Аннотации зон (КРУПНЕЕ)
plt.text(
lims[0] + 0.1, lims[1] - 0.1,
"ЗОНА СМЕРТЕЛЬНОСТИ\n(Риск смерти растет быстрее,\nчем риск травм)",
fontsize=17, # <--- Было 11, стало 15
color='#B50E3B', alpha=0.7, ha='left', va='top', fontweight='bold'
)
plt.text(
lims[1] - 0.1, lims[0] + 0.1,
"ЗОНА ТРАВМАТИЗМА\n(Высокий риск травм,\nно смерть наступает реже)",
fontsize=17, # <--- Было 11, стало 15
color='#0261C7', alpha=0.7, ha='right', va='bottom', fontweight='bold'
)
# 8. Оси и Заголовок (КРУПНЕЕ)
plt.xlabel('Влияние на переход: Легкие → Травмы', fontsize=18, fontweight='bold', color='#2d3436', labelpad=17)
plt.ylabel('Влияние на переход: Травмы → Смерть', fontsize=18, fontweight='bold', color='#2d3436', labelpad=17)
plt.title('Проверка на пропорциональность шансов', fontsize=24, fontweight='heavy', color='#2d3436', pad=25)
# Увеличиваем цифры на осях
plt.tick_params(axis='both', which='major', labelsize=16)
# Сетка и рамки
plt.grid(True, linestyle=':', color='grey', alpha=0.3)
sns.despine()
# Легенда
plt.legend(frameon=True, loc='lower right', fontsize=15)
plt.tight_layout()
plt.show()
display(comparison_df.head(10))
Обучаем бинарную модель 1 (Легкие vs Травмы/Смерть)...
/Users/skuril/.venv/lib/python3.13/site-packages/sklearn/linear_model/_linear_loss.py:200: RuntimeWarning: divide by zero encountered in matmul raw_prediction = X @ weights + intercept /Users/skuril/.venv/lib/python3.13/site-packages/sklearn/linear_model/_linear_loss.py:200: RuntimeWarning: overflow encountered in matmul raw_prediction = X @ weights + intercept /Users/skuril/.venv/lib/python3.13/site-packages/sklearn/linear_model/_linear_loss.py:200: RuntimeWarning: invalid value encountered in matmul raw_prediction = X @ weights + intercept /Users/skuril/.venv/lib/python3.13/site-packages/sklearn/linear_model/_linear_loss.py:330: RuntimeWarning: divide by zero encountered in matmul grad[:n_features] = X.T @ grad_pointwise + l2_reg_strength * weights /Users/skuril/.venv/lib/python3.13/site-packages/sklearn/linear_model/_linear_loss.py:330: RuntimeWarning: overflow encountered in matmul grad[:n_features] = X.T @ grad_pointwise + l2_reg_strength * weights /Users/skuril/.venv/lib/python3.13/site-packages/sklearn/linear_model/_linear_loss.py:330: RuntimeWarning: invalid value encountered in matmul grad[:n_features] = X.T @ grad_pointwise + l2_reg_strength * weights
Обучаем бинарную модель 2 (Выживание vs Смерть)...
/Users/skuril/.venv/lib/python3.13/site-packages/sklearn/linear_model/_linear_loss.py:200: RuntimeWarning: divide by zero encountered in matmul raw_prediction = X @ weights + intercept /Users/skuril/.venv/lib/python3.13/site-packages/sklearn/linear_model/_linear_loss.py:200: RuntimeWarning: overflow encountered in matmul raw_prediction = X @ weights + intercept /Users/skuril/.venv/lib/python3.13/site-packages/sklearn/linear_model/_linear_loss.py:200: RuntimeWarning: invalid value encountered in matmul raw_prediction = X @ weights + intercept /Users/skuril/.venv/lib/python3.13/site-packages/sklearn/linear_model/_linear_loss.py:330: RuntimeWarning: divide by zero encountered in matmul grad[:n_features] = X.T @ grad_pointwise + l2_reg_strength * weights /Users/skuril/.venv/lib/python3.13/site-packages/sklearn/linear_model/_linear_loss.py:330: RuntimeWarning: overflow encountered in matmul grad[:n_features] = X.T @ grad_pointwise + l2_reg_strength * weights /Users/skuril/.venv/lib/python3.13/site-packages/sklearn/linear_model/_linear_loss.py:330: RuntimeWarning: invalid value encountered in matmul grad[:n_features] = X.T @ grad_pointwise + l2_reg_strength * weights
| Коэфф (1 vs 2+3) | Коэфф (1+2 vs 3) | Разница | |
|---|---|---|---|
| district_Северо-Кавказский | 0.821361 | 0.319458 | 0.501903 |
| wrong_way | 0.462329 | 0.722238 | 0.259909 |
| district_Южный | 0.477841 | 0.229927 | 0.247914 |
| female_driver | -0.348171 | -0.591082 | 0.242911 |
| road_surface_cat_6 | 0.024725 | -0.209180 | 0.233905 |
| road_defects_cat_4 | -0.050776 | 0.155795 | 0.206570 |
| impaired_driving | 0.608675 | 0.796796 | 0.188121 |
| road_defects_cat_7 | 0.417215 | 0.599050 | 0.181835 |
| district_Уральский | -0.115036 | 0.065883 | 0.180919 |
| road_surface_cat_1 | -0.154668 | -0.334737 | 0.180070 |
In [52]:
# --- 1. Настройки стиля ---
try:
plt.rcParams['font.family'] = 'Montserrat'
except:
plt.rcParams['font.family'] = 'sans-serif'
plt.rcParams['font.sans-serif'] = ['Arial', 'Helvetica', 'Verdana']
sns.set_style("whitegrid")
# --- 2. Ваши переменные ---
categorical_cols = [
'district', 'road_defects_cat',
'road_surface_cat', 'lighting_cat', 'SEASON'
]
numerical_cols = [
'n_VEHICLES', 'n_guilty', 'guilty_exp_avg',
'vehicle_failure', 'female_driver', 'no_seatbelt_injury',
'impaired_driving', 'wrong_way'
]
all_vars = numerical_cols + categorical_cols
# --- 3. Вспомогательные функции ---
def get_logits(subset):
"""Считает логиты для группы"""
# P(Y <= 1)
p1 = (subset['severity'] == 1).mean()
# P(Y <= 2)
p2 = (subset['severity'] <= 2).mean()
# Защита от деления на ноль
epsilon = 1e-3
p1 = np.clip(p1, epsilon, 1 - epsilon)
p2 = np.clip(p2, epsilon, 1 - epsilon)
return np.array([np.log(p1 / (1 - p1)), np.log(p2 / (1 - p2))])
def find_intersection(y_base, y_target):
"""Находит координаты пересечения отрезков"""
yb1, yb2 = y_base
yt1, yt2 = y_target
diff1 = yt1 - yb1
diff2 = yt2 - yb2
# Если знаки разные — было пересечение
if np.sign(diff1) != np.sign(diff2) and diff1 != 0 and diff2 != 0:
m_base = yb2 - yb1
m_target = yt2 - yt1
x_cross = 1 + (yb1 - yt1) / (m_target - m_base)
y_cross = m_base * (x_cross - 1) + yb1
return x_cross, y_cross
return None
def analyze_variable(data, feature):
"""Анализирует переменную и готовит данные для графика"""
# 1. Определяем "Опасную" группу и формируем ЧИТАЕМОЕ название
if data[feature].nunique() == 2:
target_mask = data[feature] == 1
# Название: impaired_driving = 1
label_title = f"{feature} = 1"
elif feature in categorical_cols or data[feature].dtype == 'object':
# Ищем категорию с макс. смертностью
worst_cat = data.groupby(feature)['severity'].apply(lambda x: (x==3).mean()).idxmax()
target_mask = data[feature] == worst_cat
# Название: district = Южный
label_title = f"{feature} = {worst_cat}"
else:
# Численные
threshold = data[feature].quantile(0.8)
target_mask = data[feature] > threshold
# Название: High n_VEHICLES
label_title = f"High {feature}"
# 2. Считаем логиты
l_base = get_logits(data[~target_mask])
l_target = get_logits(data[target_mask])
# 3. Метрика Gap (расстояние между линиями)
gap = np.mean(np.abs(l_target - l_base))
return {
'var': feature, # Исходное имя переменной
'title': label_title, # Красивый заголовок для графика
'gap': gap,
'logits_base': l_base,
'logits_target': l_target
}
# --- 4. Расчет и Сортировка ---
results = []
for var in all_vars:
results.append(analyze_variable(df, var))
# Сортируем: Сначала большие разрывы (Group 1), в конце слипающиеся (Group 3)
sorted_results = sorted(results, key=lambda x: x['gap'], reverse=True)
# Делим на группы
group_1 = sorted_results[:5] # Топ-5 (Далекие линии)
group_2 = sorted_results[5:9] # Средние 5
group_3 = sorted_results[9:13] # Боттом-4 (Близкие/Пересекающиеся)
# --- Обновленная функция с ЕДИНЫМ МАСШТАБОМ Y ---
def plot_clean_grid_shared_axis(group_data, title, desc, n_charts):
# 1. Сначала вычисляем общие границы (min и max) для всей группы
all_values = []
for item in group_data:
all_values.extend(item['logits_base'])
all_values.extend(item['logits_target'])
y_min_global = min(all_values)
y_max_global = max(all_values)
# Добавляем отступы сверху и снизу (10% от размаха), чтобы точки не прилипали к краям
y_range = y_max_global - y_min_global
margin = y_range * 0.15
ylim_shared = (y_min_global - margin, y_max_global + margin)
# 2. Настройка сетки
n_cols = 3 if n_charts > 4 else 2
n_rows = int(np.ceil(n_charts / n_cols))
fig, axes = plt.subplots(n_rows, n_cols, figsize=(16, 6 * n_rows), dpi=120)
axes = axes.flatten()
for i, item in enumerate(group_data):
ax = axes[i]
y_base = item['logits_base']
y_target = item['logits_target']
# --- ПРИМЕНЯЕМ ЕДИНЫЙ МАСШТАБ ---
ax.set_ylim(ylim_shared)
# ЛИНИИ
ax.plot([1, 2], y_base, color='#95a5a6', linestyle='--', linewidth=2, marker='o', label='Остальные')
ax.plot([1, 2], y_target, color='#B50E3B', linestyle='-', linewidth=3.5, marker='o', label='Целевая группа')
# ПЕРЕСЕЧЕНИЕ
cross_point = find_intersection(y_base, y_target)
if cross_point:
cx, cy = cross_point
# Жирная красная точка
ax.plot(cx, cy, marker='o', markersize=14, color='#B50E3B', zorder=10)
# Надпись КАПСОМ
# Считаем отступ относительно общего масштаба, чтобы текст не прыгал
text_offset_y = (ylim_shared[1] - ylim_shared[0]) * 0.05
ax.text(cx - 0.08, cy + text_offset_y, "ПЕРЕСЕЧЕНИЕ",
color='#B50E3B', fontsize=11, fontweight='black',
ha='center', va='bottom', zorder=11)
# ОФОРМЛЕНИЕ
ax.set_title(item['title'], fontsize=14, fontweight='bold', color='#2d3436')
ax.set_xticks([1, 2])
ax.set_xticklabels(['Порог 1|2', 'Порог 2|3'], fontsize=11)
# Подпись оси Y только слева (но цифры оставляем везде для удобства)
if i % n_cols == 0:
ax.set_ylabel('Logit (Кумулятивный шанс)', fontsize=10, color='gray')
ax.grid(True, linestyle=':', alpha=0.6)
# Легенда
ax.legend(fontsize=9, loc='upper left', frameon=True)
# Удаляем пустые графики
for j in range(i + 1, len(axes)):
fig.delaxes(axes[j])
plt.suptitle(f"{title}\n{desc}", fontsize=24, fontweight='heavy', y=1.02, color='#2d3436')
plt.tight_layout()
plt.show()
# --- ЗАПУСК ПО ГРУППАМ ---
# Группа 1: Сильное влияние (Линии далеко)
plot_clean_grid_shared_axis(group_1,
"Диагностические графики пропорциональности шансов",
"",
5)
# Группа 2: Умеренное влияние
plot_clean_grid_shared_axis(group_2,
"Диагностические графики пропорциональности шансов",
"",
4)
# Группа 3: Близкие линии / Пересечения
plot_clean_grid_shared_axis(group_3,
"Диагностические графики пропорциональности шансов",
"",
4)
Обучение модели¶
In [29]:
X = X.astype(float)
model = OrderedModel(y, X, distr='logit')
res = model.fit(method='bfgs', disp=False)
print(res.summary())
OrderedModel Results
==============================================================================
Dep. Variable: severity Log-Likelihood: -4.1057e+05
Model: OrderedModel AIC: 8.212e+05
Method: Maximum Likelihood BIC: 8.217e+05
Date: Mon, 15 Dec 2025
Time: 20:47:08
No. Observations: 419576
Df Residuals: 419531
Df Model: 43
==============================================================================================
coef std err z P>|z| [0.025 0.975]
----------------------------------------------------------------------------------------------
n_VEHICLES 0.0415 0.003 13.587 0.000 0.036 0.048
n_guilty 0.1310 0.003 43.446 0.000 0.125 0.137
guilty_exp_avg 0.0216 0.003 7.118 0.000 0.016 0.027
vehicle_failure 0.2736 0.011 24.610 0.000 0.252 0.295
female_driver -0.3893 0.008 -50.681 0.000 -0.404 -0.374
no_seatbelt_injury 0.1049 0.006 17.131 0.000 0.093 0.117
impaired_driving 0.6681 0.009 74.868 0.000 0.651 0.686
wrong_way 0.5249 0.006 82.315 0.000 0.512 0.537
district_Дальневосточный -0.1242 0.014 -8.862 0.000 -0.152 -0.097
district_Приволжский -0.0639 0.009 -7.124 0.000 -0.082 -0.046
district_Северо-Западный -0.1272 0.011 -11.337 0.000 -0.149 -0.105
district_Северо-Кавказский 0.6103 0.013 48.718 0.000 0.586 0.635
district_Сибирский 0.0441 0.011 3.993 0.000 0.022 0.066
district_Уральский -0.0716 0.013 -5.646 0.000 -0.096 -0.047
district_Южный 0.3977 0.010 39.129 0.000 0.378 0.418
road_defects_cat_0 0.4604 0.114 4.039 0.000 0.237 0.684
road_defects_cat_1 0.3419 0.019 18.222 0.000 0.305 0.379
road_defects_cat_10 0.4560 0.212 2.153 0.031 0.041 0.871
road_defects_cat_11 0.2244 0.043 5.238 0.000 0.140 0.308
road_defects_cat_12 0.2874 0.029 9.792 0.000 0.230 0.345
road_defects_cat_13 0.3477 0.042 8.279 0.000 0.265 0.430
road_defects_cat_14 0.2409 0.011 21.443 0.000 0.219 0.263
road_defects_cat_2 0.2760 0.018 15.673 0.000 0.241 0.311
road_defects_cat_3 0.2757 0.318 0.866 0.386 -0.348 0.899
road_defects_cat_4 0.0095 0.094 0.101 0.919 -0.175 0.194
road_defects_cat_6 0.5877 0.028 20.875 0.000 0.533 0.643
road_defects_cat_7 0.4707 0.028 16.988 0.000 0.416 0.525
road_defects_cat_8 0.4681 0.021 22.550 0.000 0.427 0.509
road_defects_cat_9 0.4913 0.057 8.671 0.000 0.380 0.602
road_surface_cat_0 -0.0603 0.028 -2.187 0.029 -0.114 -0.006
road_surface_cat_1 -0.1983 0.015 -13.398 0.000 -0.227 -0.169
road_surface_cat_2 0.0090 0.008 1.109 0.267 -0.007 0.025
road_surface_cat_3 -0.5144 0.440 -1.169 0.242 -1.377 0.348
road_surface_cat_4 -0.0285 0.012 -2.291 0.022 -0.053 -0.004
road_surface_cat_5 -0.0917 0.070 -1.318 0.187 -0.228 0.045
road_surface_cat_6 -0.0354 0.051 -0.699 0.485 -0.135 0.064
lighting_cat_1 0.6054 0.510 1.187 0.235 -0.394 1.605
lighting_cat_2 0.3276 0.007 43.855 0.000 0.313 0.342
lighting_cat_3 -0.0411 0.011 -3.858 0.000 -0.062 -0.020
lighting_cat_4 0.1128 0.017 6.773 0.000 0.080 0.145
SEASON_1 -0.0386 0.010 -3.884 0.000 -0.058 -0.019
SEASON_2 -0.0381 0.009 -4.333 0.000 -0.055 -0.021
SEASON_4 -0.0063 0.008 -0.771 0.441 -0.022 0.010
1/2 0.2365 0.009 26.881 0.000 0.219 0.254
2/3 0.6667 0.002 294.654 0.000 0.662 0.671
==============================================================================================
Отображение результатов¶
In [30]:
num_thresholds = len(y.unique()) - 1
params = res.params[:-num_thresholds]
conf = res.conf_int().iloc[:-num_thresholds]
conf.columns = ['Lower', 'Upper']
results_df = pd.DataFrame({
'coef': params,
'lower': conf['Lower'],
'upper': conf['Upper']
})
# Считаем Odds Ratio
results_df['Odds_Ratio'] = np.exp(results_df['coef'])
# --- ГРАФИК: Топ-20 самых влиятельных факторов ---
# Сортируем по модулю коэффициента (самые сильные отклонения в любую сторону)
results_df['abs_coef'] = results_df['coef'].abs()
top_factors = results_df.sort_values(by='abs_coef', ascending=True).tail(20)
plt.figure(figsize=(12, 12))
# Рисуем
plt.errorbar(
top_factors['coef'],
top_factors.index,
xerr=[top_factors['coef'] - top_factors['lower'],
top_factors['upper'] - top_factors['coef']],
fmt='o', color='#2E8B57', ecolor='gray', capsize=3
)
plt.axvline(x=0, color='black', linestyle='--', linewidth=1)
plt.title('Топ-20 факторов, влияющих на тяжесть ДТП', fontsize=15)
plt.xlabel('← Снижает тяжесть (Безопаснее) | Повышает тяжесть (Опаснее) →', fontsize=12)
plt.grid(True, alpha=0.3)
plt.show()
# --- Вывод Топ-5 Опасных и Безопасных ---
print("\n🔥 ТОП-5 Факторов, ПОВЫШАЮЩИХ риск тяжких последствий:")
print(results_df.sort_values(by='Odds_Ratio', ascending=False)[['Odds_Ratio']].head(5))
print("\n🛡️ ТОП-5 Факторов, СНИЖАЮЩИХ риск (делают ДТП легче):")
print(results_df.sort_values(by='Odds_Ratio', ascending=True)[['Odds_Ratio']].head(5))
🔥 ТОП-5 Факторов, ПОВЫШАЮЩИХ риск тяжких последствий:
Odds_Ratio
impaired_driving 1.950499
district_Северо-Кавказский 1.840937
lighting_cat_1 1.832033
road_defects_cat_6 1.799826
wrong_way 1.690238
🛡️ ТОП-5 Факторов, СНИЖАЮЩИХ риск (делают ДТП легче):
Odds_Ratio
road_surface_cat_3 0.597858
female_driver 0.677555
road_surface_cat_1 0.820100
district_Северо-Западный 0.880530
district_Дальневосточный 0.883222
In [31]:
num_thresholds = len(y.unique()) - 1
params = res.params[:-num_thresholds]
conf = res.conf_int().iloc[:-num_thresholds]
conf.columns = ['Lower_Coef', 'Upper_Coef']
pvalues = res.pvalues[:-num_thresholds]
# 2. Собираем DataFrame
eval_df = pd.DataFrame({
'Coef': params,
'P-value': pvalues,
'Lower_Coef': conf['Lower_Coef'],
'Upper_Coef': conf['Upper_Coef']
})
# 3. Переводим в Odds Ratio (Отношение шансов)
# Это понятнее: "Риск выше в X раз"
eval_df['OR'] = np.exp(eval_df['Coef'])
eval_df['OR_Lower'] = np.exp(eval_df['Lower_Coef']) # Нижняя граница риска (в разах)
eval_df['OR_Upper'] = np.exp(eval_df['Upper_Coef']) # Верхняя граница риска (в разах)
# 4. Считаем ширину интервала (меру неопределенности)
eval_df['Uncertainty'] = eval_df['OR_Upper'] - eval_df['OR_Lower']
# 5. Категоризация значимости
def get_status(row):
if row['P-value'] > 0.05:
return '❌ Мусор (Незначимо)'
# Если интервал пересекает 1.0 (для OR), это тоже незначимо, но P-value это уже отловил
return '✅ Значимо'
eval_df['Status'] = eval_df.apply(get_status, axis=1)
# Сортируем по силе влияния (Odds Ratio)
eval_df = eval_df.sort_values(by='OR', ascending=False)
# Выводим Топ факторов с интервалами
print("Таблица влияния с учетом доверительных интервалов:")
# Показываем только значимые
display(eval_df[eval_df['Status'] == '✅ Значимо'][['OR', 'OR_Lower', 'OR_Upper', 'Uncertainty']].head(10))
Таблица влияния с учетом доверительных интервалов:
| OR | OR_Lower | OR_Upper | Uncertainty | |
|---|---|---|---|---|
| impaired_driving | 1.950499 | 1.916682 | 1.984913 | 0.068231 |
| district_Северо-Кавказский | 1.840937 | 1.796289 | 1.886694 | 0.090405 |
| road_defects_cat_6 | 1.799826 | 1.703204 | 1.901929 | 0.198724 |
| wrong_way | 1.690238 | 1.669246 | 1.711494 | 0.042249 |
| road_defects_cat_9 | 1.634395 | 1.462610 | 1.826356 | 0.363746 |
| road_defects_cat_7 | 1.601111 | 1.516480 | 1.690466 | 0.173985 |
| road_defects_cat_8 | 1.596914 | 1.533250 | 1.663221 | 0.129971 |
| road_defects_cat_0 | 1.584712 | 1.267451 | 1.981388 | 0.713937 |
| road_defects_cat_10 | 1.577737 | 1.041736 | 2.389527 | 1.347791 |
| district_Южный | 1.488351 | 1.458997 | 1.518294 | 0.059297 |
In [78]:
summary_df = pd.DataFrame({
'Коэффициент (Log-Odds)': res.params.round(3)[:-2],
'P-value': res.pvalues.round(3)[:-2],
'Odds Ratio (Шансы)': np.exp(res.params[:-2]) # Экспонента от коэффициента
})
# 2. Добавляем столбец с "человеческим" вердиктом о значимости
summary_df['Значимо?'] = summary_df['P-value'].apply(lambda x: '✅ ДА' if x < 0.05 else '❌ нет')
# 3. Сортируем по убыванию коэффициента (сверху самые опасные, снизу самые безопасные)
summary_df = summary_df.sort_values(by='Коэффициент (Log-Odds)', ascending=False)
summary_df.to_csv('sex.csv')
In [79]:
summary_df
Out[79]:
| Коэффициент (Log-Odds) | P-value | Odds Ratio (Шансы) | Значимо? | |
|---|---|---|---|---|
| impaired_driving | 0.668 | 0.000 | 1.950499 | ✅ ДА |
| district_Северо-Кавказский | 0.610 | 0.000 | 1.840937 | ✅ ДА |
| lighting_cat_1 | 0.605 | 0.235 | 1.832033 | ❌ нет |
| road_defects_cat_6 | 0.588 | 0.000 | 1.799826 | ✅ ДА |
| wrong_way | 0.525 | 0.000 | 1.690238 | ✅ ДА |
| road_defects_cat_9 | 0.491 | 0.000 | 1.634395 | ✅ ДА |
| road_defects_cat_7 | 0.471 | 0.000 | 1.601111 | ✅ ДА |
| road_defects_cat_8 | 0.468 | 0.000 | 1.596914 | ✅ ДА |
| road_defects_cat_0 | 0.460 | 0.000 | 1.584712 | ✅ ДА |
| road_defects_cat_10 | 0.456 | 0.031 | 1.577737 | ✅ ДА |
| district_Южный | 0.398 | 0.000 | 1.488351 | ✅ ДА |
| road_defects_cat_13 | 0.348 | 0.000 | 1.415849 | ✅ ДА |
| road_defects_cat_1 | 0.342 | 0.000 | 1.407677 | ✅ ДА |
| lighting_cat_2 | 0.328 | 0.000 | 1.387678 | ✅ ДА |
| road_defects_cat_12 | 0.287 | 0.000 | 1.333023 | ✅ ДА |
| road_defects_cat_2 | 0.276 | 0.000 | 1.317847 | ✅ ДА |
| road_defects_cat_3 | 0.276 | 0.386 | 1.317477 | ❌ нет |
| vehicle_failure | 0.274 | 0.000 | 1.314658 | ✅ ДА |
| road_defects_cat_14 | 0.241 | 0.000 | 1.272374 | ✅ ДА |
| road_defects_cat_11 | 0.224 | 0.000 | 1.251563 | ✅ ДА |
| n_guilty | 0.131 | 0.000 | 1.139994 | ✅ ДА |
| lighting_cat_4 | 0.113 | 0.000 | 1.119369 | ✅ ДА |
| no_seatbelt_injury | 0.105 | 0.000 | 1.110590 | ✅ ДА |
| district_Сибирский | 0.044 | 0.000 | 1.045082 | ✅ ДА |
| n_VEHICLES | 0.042 | 0.000 | 1.042395 | ✅ ДА |
| guilty_exp_avg | 0.022 | 0.000 | 1.021793 | ✅ ДА |
| road_defects_cat_4 | 0.010 | 0.919 | 1.009579 | ❌ нет |
| road_surface_cat_2 | 0.009 | 0.267 | 1.009041 | ❌ нет |
| SEASON_4 | -0.006 | 0.441 | 0.993759 | ❌ нет |
| road_surface_cat_4 | -0.029 | 0.022 | 0.971898 | ✅ ДА |
| road_surface_cat_6 | -0.035 | 0.485 | 0.965174 | ❌ нет |
| SEASON_2 | -0.038 | 0.000 | 0.962604 | ✅ ДА |
| SEASON_1 | -0.039 | 0.000 | 0.962162 | ✅ ДА |
| lighting_cat_3 | -0.041 | 0.000 | 0.959729 | ✅ ДА |
| road_surface_cat_0 | -0.060 | 0.029 | 0.941459 | ✅ ДА |
| district_Приволжский | -0.064 | 0.000 | 0.938074 | ✅ ДА |
| district_Уральский | -0.072 | 0.000 | 0.930942 | ✅ ДА |
| road_surface_cat_5 | -0.092 | 0.187 | 0.912413 | ❌ нет |
| district_Дальневосточный | -0.124 | 0.000 | 0.883222 | ✅ ДА |
| district_Северо-Западный | -0.127 | 0.000 | 0.880530 | ✅ ДА |
| road_surface_cat_1 | -0.198 | 0.000 | 0.820100 | ✅ ДА |
| female_driver | -0.389 | 0.000 | 0.677555 | ✅ ДА |
| road_surface_cat_3 | -0.514 | 0.242 | 0.597858 | ❌ нет |
Обработанные результаты¶
In [41]:
# --- 1. СЛОВАРИ И ПЕРЕВОД ---
# Свет (по порядку LabelEncoder: День, Не указано, Ночь без, Ночь с, Сумерки)
lighting_d = {
0: 'День',
1: 'Не указано',
2: 'Ночь (нет освещения)',
3: 'Ночь (с освещением)',
4: 'Сумерки'
}
# Покрытие
surface_d = {
0: 'Гололедица',
1: 'Снежный накат/Снег',
2: 'Мокрое',
3: 'Не установлено',
4: 'Обработано реагентами',
5: 'Свежий ремонт',
6: 'Грязное/Пыльное',
7: 'Сухое'
}
# Дефекты (сократил до сути, чтобы влезло на график)
road_defects_d = {
0: 'Ограничение видимости',
1: 'Плохие/Отсутствующие знаки',
2: 'Плохая зимняя уборка',
3: 'Проблемы с люками/ливневкой',
4: 'Нарушения рекламы',
5: 'Нет дефектов / Не уст.',
6: 'Плохое состояние обочин',
7: 'Проблемы с ограждениями',
8: 'Нет освещения/светофора',
9: 'Нет ТСОД в местах работ',
10: 'Дефекты Ж/Д переезда',
11: 'Нет тротуаров/остановок',
12: 'Ямы / Дефекты покрытия',
13: 'Иные недостатки',
14: 'Разметка не соотв. требованиям'
}
# Прочие переменные
static_map = {
'impaired_driving': 'Вождение в нетрезвом виде',
'wrong_way': 'Выезд на встречную',
'female_driver': 'Женщина за рулем',
'vehicle_failure': 'Техническая неисправность',
'no_seatbelt_injury': 'Непристегнутый ремень',
'n_guilty': 'Кол-во виновников',
'n_VEHICLES': 'Кол-во машин',
'guilty_exp_avg': 'Стаж водителя'
}
def get_readable_label(raw_name):
# 1. Статические имена
if raw_name in static_map:
return static_map[raw_name]
# 2. Округа
if 'district_' in raw_name:
return f"Округ: {raw_name.replace('district_', '')}"
# 3. Сезоны
if 'SEASON_' in raw_name:
seasons = {'1': 'Зима', '2': 'Весна', '3': 'Лето', '4': 'Осень'}
val = raw_name.split('_')[-1]
return f"Сезон: {seasons.get(val, val)}"
# 4. Категории с номерами
try:
val = int(raw_name.split('_')[-1])
if 'road_defects_cat' in raw_name:
return f"Дефект: {road_defects_d.get(val, str(val))}"
if 'lighting_cat' in raw_name:
return f"Свет: {lighting_d.get(val, str(val))}"
if 'road_surface_cat' in raw_name:
return f"Покрытие: {surface_d.get(val, str(val))}"
if 'road_category' in raw_name:
return f"Категория дороги: {val}"
except:
return raw_name
return raw_name
# --- 2. ПОДГОТОВКА ДАННЫХ ---
try:
plt.rcParams['font.family'] = 'Montserrat'
except:
plt.rcParams['font.family'] = 'sans-serif'
plt.rcParams['font.sans-serif'] = ['Arial', 'Helvetica', 'Verdana']
sns.set_style("white")
# Вытаскиваем данные из модели
num_thresholds = len(df['severity'].unique()) - 1
params = res.params[:-num_thresholds]
conf = res.conf_int().iloc[:-num_thresholds]
conf.columns = ['Lower_Log', 'Upper_Log']
pvalues = res.pvalues[:-num_thresholds]
eval_df = pd.DataFrame({
'OR': np.exp(params),
'OR_Lower': np.exp(conf['Lower_Log']),
'OR_Upper': np.exp(conf['Upper_Log']),
'P-value': pvalues
})
# Переименовываем индексы
eval_df.index = [get_readable_label(idx) for idx in eval_df.index]
# Фильтруем Топ-25 самых значимых
eval_df['abs_impact'] = (eval_df['OR'] - 1).abs()
plot_df = eval_df.sort_values('OR', ascending=True)
if len(plot_df) > 25:
top_risk = plot_df.tail(15)
top_safe = plot_df.head(10)
plot_df = pd.concat([top_safe, top_risk])
# --- 3. ПОСТРОЕНИЕ ГРАФИКА ---
fig, ax = plt.subplots(figsize=(13, len(plot_df) * 0.55), dpi=120)
y_range = range(len(plot_df))
for i, idx in enumerate(plot_df.index):
row = plot_df.loc[idx]
# Цвета
if row['P-value'] > 0.05:
color = '#bdc3c7' # Серый
alpha = 0.5
fontweight = 'normal'
else:
# Малиновый (Риск) / Синий (Защита)
color = '#B50E3B' if row['OR'] > 1 else '#0261C7'
alpha = 1.0
fontweight = 'bold'
# Зебра
if i % 2 == 0:
ax.axhspan(i - 0.5, i + 0.5, color='#f7f9fa', zorder=0)
# Усы
ax.hlines(i, row['OR_Lower'], row['OR_Upper'], color=color, linewidth=2, alpha=alpha, zorder=2)
# Точка
ax.scatter(row['OR'], i, color=color, s=90, alpha=alpha, zorder=3, edgecolor='white', linewidth=0.8)
# Текст значения
ax.text(
x=plot_df['OR_Upper'].max() * 1.05,
y=i,
s=f"{row['OR']:.2f}x",
va='center', ha='left',
fontsize=11, fontweight=fontweight, color=color
)
# Вертикальная линия (База)
ax.axvline(x=1, color='#2d3436', linestyle='--', linewidth=1.5, alpha=0.8, zorder=1)
# Оси
ax.set_yticks(y_range)
ax.set_yticklabels(plot_df.index, fontsize=12, fontweight='medium', color='#2d3436')
ax.set_xlabel('Отношение шансов (Odds Ratio)\n← Снижает риск | 1.0 (База) | Повышает риск →',
fontsize=13, fontweight='bold', color='#2d3436', labelpad=15)
ax.set_title('Сила влияния факторов (доверительный интервал 95%)', fontsize=20, fontweight='heavy', color='#2d3436', pad=25)
# Чистка
sns.despine(left=True, bottom=False, right=True, top=True)
ax.tick_params(axis='y', length=0)
ax.grid(axis='x', linestyle=':', alpha=0.5)
# Границы графика (с запасом под текст справа)
plt.xlim(0, plot_df['OR_Upper'].max() * 1.25)
plt.tight_layout()
plt.show()
In [115]:
# --- 1. СЛОВАРИ И ПЕРЕВОД ---
# Свет (по порядку LabelEncoder: День, Не указано, Ночь без, Ночь с, Сумерки)
lighting_d = {
0: 'День',
1: 'Не указано',
2: 'Ночь (нет освещения)',
3: 'Ночь (с освещением)',
4: 'Сумерки'
}
# Покрытие
surface_d = {
0: 'Гололедица',
1: 'Снежный накат/Снег',
2: 'Мокрое',
3: 'Не установлено',
4: 'Обработано реагентами',
5: 'Свежий ремонт',
6: 'Грязное/Пыльное',
7: 'Сухое'
}
# Дефекты (сократил до сути, чтобы влезло на график)
road_defects_d = {
0: 'Ограничение видимости',
1: 'Плохие/Отсутствующие знаки',
2: 'Недостатки зимнего содержания',
3: 'Проблемы с люками/ливневкой',
4: 'Нарушения рекламы',
5: 'Нет дефектов / Не уст.',
6: 'Плохое состояние обочин',
7: 'Проблемы с ограждениями',
8: 'Нет освещения/светофора',
9: 'Нет ТСОД в местах работ',
10: 'Дефекты Ж/Д переезда',
11: 'Нет тротуаров/остановок',
12: 'Ямы / Дефекты покрытия',
13: 'Иные недостатки',
14: 'Разметка не соотв. требованиям'
}
# Прочие переменные
static_map = {
'impaired_driving': 'Вождение в нетрезвом виде',
'wrong_way': 'Выезд на встречную',
'female_driver': 'Женщина за рулем',
'vehicle_failure': 'Техническая неисправность',
'no_seatbelt_injury': 'Непристегнутый ремень',
'n_guilty': 'Кол-во виновников',
'n_VEHICLES': 'Кол-во машин',
'guilty_exp_avg': 'Стаж водителя'
}
def get_readable_label(raw_name):
# 1. Статические имена
if raw_name in static_map:
return static_map[raw_name]
# 2. Округа
if 'district_' in raw_name:
return f"Округ: {raw_name.replace('district_', '')}"
# 3. Сезоны
if 'SEASON_' in raw_name:
seasons = {'1': 'Зима', '2': 'Весна', '3': 'Лето', '4': 'Осень'}
val = raw_name.split('_')[-1]
return f"Сезон: {seasons.get(val, val)}"
# 4. Категории с номерами
try:
val = int(raw_name.split('_')[-1])
if 'road_defects_cat' in raw_name:
return f"Дефект: {road_defects_d.get(val, str(val))}"
if 'lighting_cat' in raw_name:
return f"Свет: {lighting_d.get(val, str(val))}"
if 'road_surface_cat' in raw_name:
return f"Покрытие: {surface_d.get(val, str(val))}"
if 'road_category' in raw_name:
return f"Категория дороги: {val}"
except:
return raw_name
return raw_name
# --- 2. ПОДГОТОВКА ДАННЫХ ---
try:
plt.rcParams['font.family'] = 'Montserrat'
except:
plt.rcParams['font.family'] = 'sans-serif'
plt.rcParams['font.sans-serif'] = ['Arial', 'Helvetica', 'Verdana']
sns.set_style("white")
# --- 0) (ваш код подготовки словарей/labeling/обучения модели — без изменений) ---
# --- 1) Вытаскиваем данные из модели (ваш блок — без изменений) ---
num_thresholds = len(df['severity'].unique()) - 1
params = res.params[:-num_thresholds]
conf = res.conf_int().iloc[:-num_thresholds]
conf.columns = ['Lower_Log', 'Upper_Log']
pvalues = res.pvalues[:-num_thresholds]
eval_df = pd.DataFrame({
'OR': np.exp(params),
'OR_Lower': np.exp(conf['Lower_Log']),
'OR_Upper': np.exp(conf['Upper_Log']),
'P-value': pvalues
})
# Переименовываем индексы (ваша функция get_readable_label)
eval_df.index = [get_readable_label(idx) for idx in eval_df.index]
# --- 2) Функция отрисовки одного forest-plot ---
def plot_forest(plot_df, title, x_max=None):
plot_df = plot_df.sort_values('OR', ascending=True)
try:
plt.rcParams['font.family'] = 'Montserrat'
except:
plt.rcParams['font.family'] = 'sans-serif'
plt.rcParams['font.sans-serif'] = ['Arial', 'Helvetica', 'Verdana']
sns.set_style("white")
fig, ax = plt.subplots(figsize=(13, max(4, len(plot_df) * 0.55)), dpi=120)
y_range = range(len(plot_df))
for i, idx in enumerate(plot_df.index):
row = plot_df.loc[idx]
# Цвета
if row['P-value'] > 0.05:
color = '#bdc3c7' # Серый
alpha = 0.5
fontweight = 'normal'
else:
color = '#B50E3B' if row['OR'] > 1 else '#0261C7' # Риск / Защита
alpha = 1.0
fontweight = 'bold'
# Зебра
if i % 2 == 0:
ax.axhspan(i - 0.5, i + 0.5, color='#f7f9fa', zorder=0)
# Усы (CI)
ax.hlines(i, row['OR_Lower'], row['OR_Upper'], color=color, linewidth=2, alpha=alpha, zorder=2)
# Точка (OR)
ax.scatter(row['OR'], i, color=color, s=90, alpha=alpha, zorder=3,
edgecolor='white', linewidth=0.8)
# Текст значения
right_text_x = (x_max if x_max is not None else plot_df['OR_Upper'].max()) * 1.05
ax.text(
x=right_text_x,
y=i,
s=f"{row['OR']:.2f}x",
va='center', ha='left',
fontsize=11, fontweight=fontweight, color=color
)
# Вертикальная линия (База)
ax.axvline(x=1, color='#2d3436', linestyle='--', linewidth=1.5, alpha=0.8, zorder=1)
# Оси/подписи
ax.set_yticks(list(y_range))
ax.set_yticklabels(plot_df.index, fontsize=12, fontweight='medium', color='#2d3436')
ax.set_xlabel(
'Отношение шансов (Odds Ratio)\n← Снижает риск | 1.0 (База) | Повышает риск →',
fontsize=13, fontweight='bold', color='#2d3436', labelpad=15
)
ax.set_title(title, fontsize=20, fontweight='heavy', color='#2d3436', pad=25)
sns.despine(left=True, bottom=False, right=True, top=True)
ax.tick_params(axis='y', length=0)
ax.grid(axis='x', linestyle=':', alpha=0.5)
# Границы (с запасом под текст справа)
xmax = x_max if x_max is not None else plot_df['OR_Upper'].max()
ax.set_xlim(0, xmax * 1.25)
plt.tight_layout()
plt.show()
is_defect = pd.Index(eval_df.index).astype(str).str.startswith('Дефект:')
plot_df_defects = eval_df.loc[is_defect].copy()
plot_df_other = eval_df.loc[~is_defect].copy()
tmp = plot_df_other.copy()
# --- Топ‑7 и бот‑7 ТОЛЬКО для прочих факторов ---
plot_df_other = plot_df_other.sort_values('OR', ascending=True) # порядок по OR [web:35]
if len(plot_df_other) > 14:
plot_df_other = pd.concat(
[plot_df_other.head(7), plot_df_other.tail(7)]
) # берём 7 сверху и 7 снизу [web:22][web:42]
plot_df_other = plot_df_other.sort_values('OR', ascending=True) # вернуть правильный порядок [web:35]
# дефекты показываем все как есть (либо там можете своё ограничение сделать)
plot_df_defects = plot_df_defects.sort_values('OR', ascending=True)
plot_df_another = pd.concat([tmp, plot_df_other]).drop_duplicates(keep=False)
common_xmax = eval_df['OR_Upper'].max()
plot_forest(
plot_df_defects,
'Дефекты: сила влияния (доверительный интервал 95%)',
x_max=common_xmax
)
plot_forest(
plot_df_other,
'Некоторые прочие факторы: сила влияния (доверительный интервал 95%)',
x_max=common_xmax
)
plot_forest(
plot_df_another,
'Оставшиеся прочие факторы: сила влияния (доверительный интервал 95%)',
x_max=common_xmax
)
Переход к вероятностям¶
In [93]:
BASE_DEFECT_COL = "road_defects_cat_5" # так, как в ноутбуке [file:57]
BASE_DEFECT_COL_NUM = 5
road_defects_d = {
0: 'Ограничение видимости',
1: 'Плохие/Отсутствующие знаки',
2: 'Недостатки зимнего содержания',
3: 'Проблемы с люками/ливневкой',
4: 'Нарушения рекламы',
5: 'Нет дефектов / Не уст.',
6: 'Плохое состояние обочин',
7: 'Проблемы с ограждениями',
8: 'Нет освещения/светофора',
9: 'Нет ТСОД в местах работ',
10: 'Дефекты Ж/Д переезда',
11: 'Нет тротуаров/остановок',
12: 'Ямы / Дефекты покрытия',
13: 'Иные недостатки',
14: 'Разметка не соотв. требованиям'
}
# палитра (как в ноутбуке) [file:57]
hexcolors = ["#B50E3B", "#850AD6", "#490FD2", "#153AE0", "#0261C7"]
severity_palette = {1: hexcolors[-1], 2: hexcolors[2], 3: hexcolors[0]}
# --- 1. Настройки стиля (крупный шрифт, узкие графики) ---
plt.rcParams.update({
"font.family": "sans-serif",
"font.sans-serif": ["Arial", "Helvetica", "Verdana"],
"font.size": 16,
"axes.titlesize": 24,
"axes.labelsize": 18,
"xtick.labelsize": 16,
"ytick.labelsize": 16,
})
# Параметры цветов
ALPHA = 0.05
gray_color = "#BDC3C7" # Серый для незначимых
# Цвета для severity: 1 (синий), 2 (фиолетовый), 3 (красный)
hexcolors = ["#B50E3B", "#850AD6", "#490FD2", "#153AE0", "#0261C7"]
severity_palette = {1: hexcolors[-1], 2: hexcolors[2], 3: hexcolors[0]}
# --- 2. Расчет Marginal Effects (MEM) ---
# Ищем колонки дефектов (имена должны совпадать с X.columns)
defect_cols = sorted([c for c in X.columns if "roaddefectscat" in c or "road_defects_cat" in c])
if not defect_cols:
raise ValueError("Колонки дефектов не найдены! Проверьте X.columns.")
# Базовая точка: средние значения, дефекты = 0
x0 = X.mean(axis=0).copy()
x0.loc[defect_cols] = 0.0
# Функция предсказания вероятностей
def predict_probs_ordered(model_res, exog_row):
exog_df = pd.DataFrame([exog_row.values], columns=exog_row.index)
probs = model_res.model.predict(model_res.params, exog=exog_df)
return np.asarray(probs).flatten()
# Базовые вероятности (p0)
p0 = predict_probs_ordered(res, x0)
K = len(p0) # Количество классов (обычно 3)
rows = []
# Достаем p-values прямо из модели
pvalues_series = res.pvalues
for col in defect_cols:
# Симулируем наличие дефекта (1.0)
xi = x0.copy()
xi[col] = 1.0
pi = predict_probs_ordered(res, xi)
# P-value для этой переменной
pval = pvalues_series.get(col, np.nan)
for k in range(K):
rows.append({
"defect_col": col,
"severity": k + 1,
"delta_pp": 100.0 * (pi[k] - p0[k]), # Разница в % пунктах
"p_value": float(pval)
})
effects = pd.DataFrame(rows)
# --- 3. Перевод названий и сортировка ---
def get_defect_label(col_name):
import re
# Ищем цифру в конце названия колонки
m = re.search(r"(\d+)$", col_name)
if m:
num = int(m.group(1))
# road_defects_d - словарь из твоего ноутбука
if 'road_defects_d' in globals() and num in road_defects_d:
return road_defects_d[num]
return col_name
effects["defect_ru"] = effects["defect_col"].apply(get_defect_label)
# Сортировка: определяем порядок по влиянию на самый тяжелый класс (например, 3)
target_sev = K
order_df = effects[effects["severity"] == target_sev].sort_values("delta_pp", ascending=True)
order_list = order_df["defect_ru"].tolist()
# --- 4. Построение 3-х отдельных графиков ---
# --- 4. Построение 3-х отдельных графиков (с исправленным отступом) ---
for k in range(1, K + 1):
fig, ax = plt.subplots(figsize=(12, 8), dpi=150) # Сделали чуть шире (12)
subset = effects[effects["severity"] == k].copy()
subset["defect_ru"] = pd.Categorical(subset["defect_ru"], categories=order_list, ordered=True)
subset = subset.sort_values("defect_ru")
# Цвета
bar_colors = []
for _, row in subset.iterrows():
if row["p_value"] > ALPHA or pd.isna(row["p_value"]):
bar_colors.append(gray_color)
else:
bar_colors.append(severity_palette.get(k, "#333333"))
bars = ax.barh(subset["defect_ru"], subset["delta_pp"], color=bar_colors, edgecolor="none", height=0.7)
# Расширяем границы оси X, чтобы влез текст
# Находим мин и макс значения
min_val = subset["delta_pp"].min()
max_val = subset["delta_pp"].max()
# Добавляем 15-20% запаса по краям для подписей
range_span = max_val - min_val
if range_span == 0: range_span = 1
ax.set_xlim(min_val - range_span * 0.2, max_val + range_span * 0.2)
# Заголовки
severity_names = {1: "Тяжесть 1 (Легкие)", 2: "Тяжесть 2 (Тяжелые)", 3: "Тяжесть 3 (Смертельные)"}
title_text = severity_names.get(k, f"Тяжесть {k}")
ax.set_title(f"Влияние дефектов: {title_text}", pad=20, fontweight="bold")
ax.set_xlabel("Изменение вероятности (п.п.)", color="#555555", labelpad=10)
ax.axvline(0, color="black", linewidth=1, linestyle="-", alpha=0.8)
ax.grid(axis="x", linestyle="--", alpha=0.4)
sns.despine(ax=ax, left=True, bottom=False)
# Подписи
for bar, val in zip(bars, subset["delta_pp"]):
# Смещаем текст сильнее
offset = 0.5 if val >= 0 else -0.5
ha_align = 'left' if val >= 0 else 'right'
ax.text(val + offset,
bar.get_y() + bar.get_height()/2,
f"{val:+.1f}",
va='center',
ha=ha_align,
fontsize=14,
fontweight='bold',
color="black")
plt.tight_layout()
plt.show()